Visualizing the Impact of Palm Oil Production
Wanda Wang
April 24, 2016
Introduction
Deforestation is accelerating at an alarming rate around the world as a result of Palm oil production. In order to satisfy growing global demand for Palm oil, an ingredient commonly found in cosmetics to processed foods - clearing the way for lucrative oil plantations by forest fires is negatively impacting our environment. Concerns include the loss of rainforest habitat for endangered species, harmful carbon dioxide emissions, coupled with the displacement of local indigenous communities.
Exploratory Questions
What countries are the top producers of Palm oil? Where are the palm oil mills located? Where are the forest-fire hotspots?
The Datasets
The UN Food and Agriculture organization provides information on global Palm oil production and area harvested. The Global Forest Watch also has Palm oil Mill location data and fire alert data available for our review.
Industry Overview
To navigate to where Palm Oil production is at its highest volume, we first approach the UN Food and Agriculture dataset:
setwd("/Users/drscholls303/Desktop/Project1")
# Load the required R packages and specified datasets
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 3.2.4faostat.palm <-read.csv("Faostats.csv")
faostat.palm <- faostat.palm %>% na.omit()
#Select relevant columns
producers <-select(faostat.palm, AreaName, ElementName, Value, Year, ItemName)
# Filter for top producers in most recent available year of 2014
producers.year<-filter(producers, Year == 2014, ItemName == "Oil, palm") %>% group_by(AreaName, ElementName='Production') %>%
summarise(., Production=sum(Value))
nrow(producers.year) #15 rows ## [1] 15#Calculate percentage of world production for each country
producers.year$Percent = round(producers.year$Production[1:15]*100/
sum(producers.year$Production), digits = 1)
producers.year<-tbl_df(producers.year)
top<-top_n(producers.year,5,Production)
top## Source: local data frame [5 x 4]
##
## AreaName ElementName Production Percent
## (fctr) (chr) (int) (dbl)
## 1 Colombia Production 1109586 2.1
## 2 Germany Production 605582 1.1
## 3 Indonesia Production 29344500 55.6
## 4 Malaysia Production 19667016 37.2
## 5 Papua New Guinea Production 520000 1.0top$Production = top$Production/1000
top## Source: local data frame [5 x 4]
##
## AreaName ElementName Production Percent
## (fctr) (chr) (dbl) (dbl)
## 1 Colombia Production 1109.586 2.1
## 2 Germany Production 605.582 1.1
## 3 Indonesia Production 29344.500 55.6
## 4 Malaysia Production 19667.016 37.2
## 5 Papua New Guinea Production 520.000 1.0#top 5 graph
library(scales)
top.plot<-ggplot(top, aes(AreaName, Production, fill=factor(AreaName), label=top$Percent)) + geom_bar(stat="identity") + xlab("Country") + ylab("Tonnes(1000's)") + ggtitle("Palm Oil Production(%)")
top.plot<-top.plot + geom_text()+ guides(fill = guide_legend(title = "Top Five", title.position = "top"))
top.plot 
The top five producers in the world include Indonesia and Malaysia. Indonesia accounts for over 55% or worldwide production. Malaysia accounts for around 37%. The combined 92% allocation in this region should be noted.
Industry Map - Bird’s eye view
library(maptools)## Loading required package: sp## Checking rgeos availability: TRUElibrary(ggplot2)
library(ggmap)
library(rgeos)## rgeos version: 0.3-11, (SVN revision 479)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.1-0
## Polygon checking: TRUElibrary(dplyr)
plotmap <- readShapePoly(fn="TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp")
palmHa <-read.csv("Faostats.csv", stringsAsFactors=F)
palmHa <-filter(palmHa, ElementName == 'Area harvested')
countriez <-distinct(palmHa, AreaName)
countriez$AreaName[1:17]## [1] "Indonesia" "Malaysia"
## [3] "Nigeria" "Thailand"
## [5] "Ghana" "Guinea"
## [7] "Democratic Republic of the Congo" "CÌ«te d'Ivoire"
## [9] "Colombia" "Ecuador"
## [11] "Papua New Guinea" "Cameroon"
## [13] "Honduras" "Brazil"
## [15] "Costa Rica" "Guatemala"
## [17] "Philippines"countriez$AreaName[19:44]## [1] "China"
## [2] "Mexico"
## [3] "Peru"
## [4] "Venezuela (Bolivarian Republic of)"
## [5] "Benin"
## [6] "Sierra Leone"
## [7] "Paraguay"
## [8] "Angola"
## [9] "Guinea-Bissau"
## [10] "Liberia"
## [11] "Togo"
## [12] "Dominican Republic"
## [13] "Solomon Islands"
## [14] "Congo"
## [15] "Senegal"
## [16] "Central African Republic"
## [17] "Burundi"
## [18] "Panama"
## [19] "Suriname"
## [20] "United Republic of Tanzania"
## [21] "Nicaragua"
## [22] "Gabon"
## [23] "Gambia"
## [24] "Equatorial Guinea"
## [25] "Madagascar"
## [26] "Sao Tome and Principe"names(plotmap)[5] <- "AreaName"
target <- c(countriez$AreaName[1:17],countriez$AreaName[19:44])
#Use foritfy to transform spatial plotmap into dataframe
plotmapDf <- fortify(plotmap, region = "AreaName") %>%
filter(id %in% target)
palmHaMapDf <- merge(plotmapDf, palmHa, by.x="id", by.y="AreaName")
backgroundMap <- ggplot(data=palmHaMapDf) + geom_path(aes(x=long, y=lat, group=group), color='grey') + coord_equal() + geom_polygon(aes(x=long, y=lat, group=group, fill=Value))
mapAsia <- get_map(location = 'Indonesia', zoom=4)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Indonesia&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Indonesia&sensor=falsemapAfrica <- get_map(location = 'Africa', zoom=4)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Africa&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Africa&sensor=falsemapSouthAm <- get_map(location = 'South America', zoom=4)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=South+America&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=South%20America&sensor=falseggmapObjAsia <- ggmap(mapAsia)
Asia <- ggmapObjAsia + geom_polygon(aes(x=long, y=lat, group=group, fill=Value), data=palmHaMapDf, alpha=.9) +
geom_path(aes(x=long, y=lat, group=group), data=palmHaMapDf, color='black')
Asia## Warning: Removed 331938 rows containing missing values (geom_path).
ggmapObjAf <- ggmap(mapAfrica)
Africa <- ggmapObjAf + geom_polygon(aes(x=long, y=lat, group=group, fill=Value), data=palmHaMapDf, alpha=.9) +
geom_path(aes(x=long, y=lat, group=group), data=palmHaMapDf, color='black')
Africa## Warning: Removed 22194 rows containing missing values (geom_path).
ggmapSouthAm <- ggmap(mapSouthAm)
SouthAm <- ggmapSouthAm + geom_polygon(aes(x=long, y=lat, group=group, fill=Value), data=palmHaMapDf, alpha=.9) +
geom_path(aes(x=long, y=lat, group=group), data=palmHaMapDf, color='black')
SouthAm## Warning: Removed 53298 rows containing missing values (geom_path).
As a visual confirmation on the map - looking at SouthEast Asia vs Africa or South America, there are more acres harvested within Indonesia.
Palm Oil Harvested in Indonesia
library(dplyr)
library(ggplot2)
faostat.palm <-read.csv("Faostats.csv")
faostat.palm <- faostat.palm %>% na.omit()
AH<-filter(faostat.palm, ElementName=='Area harvested', AreaName=='Indonesia')
AH$Value = AH$Value/1000
areaGraph<-ggplot(AH, aes(Year, Value, fill=AreaName)) + geom_bar(stat="identity") +
ylab("Hectares(1000's)") + ggtitle("Area Harvested through Time")
areaGraph
The rise in Area Harvested over time is alarming over time. The drastic percent change from 2000 to 2014 rose from a value of 2,014,000 to 7,407,090 hectares impacted or 2.7 times.
Mill Distribution and Fire Hotspots
library(dplyr)
library(leaflet)
mills<-read.csv("mills.csv")
fires<-read.csv("fires.csv")
last.year.fires <-filter(fires, OBJECTID > 11556)#fire-alerts, as of 1/01/2015
#zoom in on SE asia, Circles for mills, Markers for fire alerts
map <- leaflet(mills) %>%
addTiles('http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png',attribution='Map tiles by <a href="http://stamen.com">Stamen Design</a>, <a href="http://creativecommons.org/licenses/by/3.0">CC BY 3.0</a> — Map data © <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>') %>% addCircles(~longitude, ~latitude, popup=mills$type, weight = 3, radius=40, color="#ffa500", stroke = TRUE, fillOpacity = 0.8) %>%
addMarkers(data = last.year.fires, ~Latitude, ~Longitude, clusterOptions = markerClusterOptions()) %>%
addLegend("bottomright", colors= "#ffa500", labels="Mills", title="Legend")
mapConclusion
There is a strong presence of fire hotspots within the last year where the majority of Mills are located. Indonesia possesses the most Palm Oil Mills as currently reported in available data.
Future Questions
Are there any geographic overlaps with existing nature reservations- what is the tree cover loss over time in those areas? How much has the Air Quality Index been impacted? Is there a trend between species extinction and habitat location in Palm Oil producing nations?